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Abstract 

tin This paper introduces a new shape-based image reconstruction technique apph- 

OO cable to a large class of imaging problems formulated in a variational sense. Given 

CN a collection of shape priors (a shape dictionary), we define our problem as choosing 

the right elements and geometrically composing them through basic set operations to 
characterize desired regions in the image. This combinatorial problem can be relaxed 
and then solved using classical descent methods. The main component of this relax- 
ation is forming certain compactly supported functions which we call "knolls", and 
reformulating the shape representation as a basis expansion in terms of such func- 
tions. To select suitable elements of the dictionary, our problem ultimately reduces 
^ to solving a nonlinear program with sparsity constraints. We provide a new sparse 

I— ' nonlinear reconstruction technique to approach this problem. The performance of 

^ proposed technique is demonstrated with some standard imaging problems including 

^ image segmentation. X-ray tomography and diffusive tomography. 

OO 
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^ 1 Introduction 

O 

CO In many imaging applications, the main objective is to identify and characterize regions 

of interest in a given domain. This characterization is usually based on some property 



defined over the imaging domain. For instance in image segmentation 40,42], this 



^ property is directly or statistically related to the pixel values and the partitioning is 



^ usually meant to aggregate similar regions. In shape-based inverse problems [18||24||32] , 



the property of interest corresponds to a spatial physical parameter which needs to be 
determined by some indirect observations. Here a shape-based characterization delineates 
inclusions and obstacles causing a contrast in the values of the spatial parameter. In 
the basic binary case the shape characterization problem is formulated as partitioning 
a compact imaging domain g M^, into (l and D \ (l^ where g is a closed set 
representing the object (s) of interest. 
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In this paper we propose a new approach to shape-based modehng. The idea is to 
characterize the target inclusion (or partition) as a composition of given shape prototypes. 
These prototypes may either be shapes with simple geometries or shape priors that are 
likely to be related to the structure of the target geometry. The idea of composing basic 
shapes to form more complex structures is quite intuitive. The challenge, is to formulate 
shape composition concretely and then incorporate it into the reconstruction process in 
a computationally tractable manner. 

For the purpose of shape composition, we consider a mechanism to merge shapes and 
exclude undesirable regions from the aggregate to generate fi. More formally, given a 
set of prototype shapes, Si,S2,'"Sni which are known closed sets in R^, we formulate 
the shape composition process as applying set union and relative complement among a 
selected number of shapes in the reference set to generate a region that approximates 
(l. A brute force approach would entail exploring all possible selections and composition 
possibilities to find a good fit. To have an idea about the accuracy of an estimate of 
Q, we specifically focus on variational shape-based problems where a score (value of an 
energy functional) is assigned to every closed region ft e D and ft is meant to minimize 
such functional. 

This type of brute force search is in general computationally intractable. In this 
paper we will propose a relaxation formed by defining so called knoll functions which 
characterize Si^S2r"^n by their support. We will show that for a smooth energy func- 
tional associated with the shape-based problem, this relaxation can convert the problem 
into minimizing a smooth nonlinear function in R^. 

To encourage our target shape to be as simple as possible, we impose sparsity con- 
straints on the resulting nonlinear cost. Accordingly, we propose a minimization tech- 
nique inspired by the idea developed by van den Berg and Friedlander in |47 . The 
original idea in |47] developed for linear inverse problems is merely applicable to the 
corresponding quadratic least squares cost. What we will put forth is a generalization of 
this approach applied to quadratic estimates of the nonlinear problem at iterative stages 
of approaching a minima. 

The organization of this paper is as follows. In the remainder of the introduction 
section we provide a rather extensive overview of variational shape-based methods, spe- 
cially in context of active contours, pixel based and parametric level set techniques. After 
providing this background, in Section 2 we present the shape composition idea and dis- 
cuss employing a so called pseudo-logical shape interaction property which will assist 
us in developing a relaxation to the corresponding combinatorial problem. Section 3 is 
devoted to developing a Gauss-Newton type sparsity promoting algorithm to solve the 
resulting nonlinear problem. Finally, in Section 4 we consider applying the proposed 
technique to some image segmentation and shape-based inverse problem^ The image 



^The Matlab code for the proposed algorithm is available at http : //users . ece . gatech . edu/ 



processing applications considered are segmentation with missing pixels and machine| 
text recognition. We also consider an example of medical X-ray computed tomo graphy] 
[and an archaeologic resistance tomography problem, both with limited available data,] 
[whe re majo rity of the state of the art techniques would not be able to provide satisfying! 
Irecon struc tions. I 



1.1 Background on Variational Shape Reconstruction 



In the past decades considerable effort has been devoted to explorin 


g variational tech- 


niques applicable to the shape-based characterization problem |15, 


16, 


23,31 


,36146]. The 


main intuition behind a variational approach is forming an energy functional, minimiz- 


ing which tends to solve the characterization problem. More specifically, for 


an arbitrary 


bounded set g D, 




= argmin 






(1) 



where £ is the underlying energy or cost functional. An example of £ in the context of 



inverse problems is the data-model mismatch functional 



mn) = \\v - M(n)h7, 



where v represents the observed data, Ai denotes the model mapping the shape geometry! 
to the observations and S>y is a Hilbert space associated with the data. One of the main 



motivations in representing problems in variational forms is the chance of using descent 



optimization techniques. 



Of course a straight-forward approach in determining a region is determining its 
boundary C. An early technique of this type is the active contour model (snakes), where 



the shape determination amounts to determinin g the par ameters as sociated with a spline 
representation of C |31|. Starting with an initial shape representation, this process takes 



steps along the descent direction of the energy to evolve the shape towards an optimum 
state. More specifically, by defining an artificial time t, an initial contour C(0) is evolved 



in time and according to 



at 



= - 'E\C)[C(t)l 



(3) 



to find the boundary of the shape that locally minimizes £ [14|. Here 'E'{C)[.] IS a 
linear operator representing the Gateaux derivative (or first order variations) of the 
energy functional with respect to C. 



Another well known shape-based technique is the level set method 



37 



39 



. The main 



advantage of level sets over earlier methods such as snakes is their topological flexibility, 
dispelling the n eed to a ny prior assumptions about the number of connected components 
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in ft. Here the zero level set of a Lipschitz continuous function, 0, is used to identify C. 
More specifically the objective is to determine such that 



<0 xeD\n ' ^ ^ 

By using a map as Q our geometric problem is cast as the calculus of variations problem 

4) = argmin£(Jl^), (5) 

where fi^ is the shape resulted at the zero level set of 0, and O represents a certain func- 
tion space to prevent ill-conditioning. Most level set implementations consider elements 
of 6 to be signed distance functions (SDFs) [37] . 

For the level set methods, minimization of ^ is performed by evolving an initial 
level set function, (j){x^t) = (j)o{x), through the Hamilton- Jacobi equation 

^^^ + V{x,t)-V<Pix,t) = 0. (6) 

This equation results from a straightforward differentiation of the front equation t) - 
with respect to t. In this equation V{x^t) - dx/dt is a speed function applied to the 
zero level set of (j){x^t) and taken in the descent direction of £ to reduce it as (j) evolves. 
At every iteration, through a proper speed function extension or re-initialization, the 



level set function is assured to remain an SDF. We refer the reader to 10,23,37] for more 
details about implementation of this technique. 

Level set methods are among the most successful techniques in shape representations, 
mainly due to their topological flexibility. They perform remarkable for a large class of 



image processing applications 20,38], however, in general there are implementation con- 
cerns and performance deficiencies associated with them. Numerically speaking, evolving 
the level set function requires discretizing it over a dense grid of pixels and updating the 
pixel values over discrete time frames (iterations) |37|. Using pixels to parameterize the 
level set function brings a large dimensionality to the problem which may severely affect 
the overall performance of the method in dealing with ill-posed functionals and inverse 
problems (see examples in fl]). Although in principle level sets are topologically flexible, 
it is usually not possible to create holes or shapes away from the boundaries of the evolv- 
ing shape |9 . Moreover, implementation complexities such as speed function extension 
and maintaining (j) as an SDF are usually the inevitable components of this technique. 

To overcome these problems and still take advantage of the topological flexibility, 
researchers have started considering parametric forms for the level set function [l]|5]. 
This strategy not only reduces problem's dimensionality by making it parametric, also 
provides a predetermined functional form to dispel the need to re-initialization. An 



example of such effort is the work by Bernard et al |5 , where according to ([5]), O is 
the space spanned by pre-assigned B-sphne functions. The problem in this case usually 
amounts to the energy minimization 

a = arg min ) , (7) 

where cx is the vector of parameters associated with the parametric level set function 
0(x,a). 

A main advantage of parameterizing the problem as ^ is casting the shape charac- 
terization problem as a classic finite dimensional minimization problem, where gradient 
descent or even Newton type methods may be applied. For sufficiently smooth func- 
tional, gradient and Hessian of £ with respect to elements of a. may be conveniently 
calculated by following the chain rule in Frechet spaces [l]. More specifically 

|^ = E'W)[|^], (8) 

and ^ ^ 

where the linear and bilinear operators and .] are respectively the first 

and second Gateaux derivatives of £ with respect to (p. We refer the reader to |1| for 
more details on a parametric level set representation. 

Although a parametric form solves many concerns with traditional level sets, the 
problem still remains on choosing a suitable parametric form (e.g., a suitable set of basis 
functions) and efficiently determining the number of terms. Clearly a dense basis set to 
increase the resolution of the reconstructions may bring redundancy and ill-conditioning 
to the problem as the case with pixel based level sets. 

In the sequel, we propose a new way of approaching the shape problem. We extend 
one of the basic ideas developed by Aghasi et al. in 1 1 and provide a new shape recon- 
struction technique applicable to a large class of problems. The method ultimately finds 
a parametric level set form, however, we link it to a systematic algorithm of choosing 
the right parameters in the course of reconstruction. 



2 Basic Shape Composition Idea 

In solving ([T]) for the optimum shape fi, an intuitive idea would be to consider a collection 
of fixed shapes and reconstruct Q by applying basic set operations on them. More 
specifically, consider the shape collection 2) = ^2, •••^nc^}? where for i = l,2,---,n(i, 
elements Si represent closed known regions (shapes) in V. We name T) as the shape 



dictionary or simply the dictionary. Suppose through an oracle we know what the true 
shape is. To logically (in the sense of set operations) express in terms of elements of 
2), one way of approaching the problem is to select an appropriate set of shapes indexed 
by Xe ^ {1,2,'", rid} apply the union operation over them to form the more bulky 
superset Uiex® foi" ^- We now start to choose shapes indexed by Iq c {l,2,---,n^} to 
carve out portions of UieX® make it a better approximation to (in the sense 

of reducing £). To be concise, for a given dictionary D we define the objective as the 
combinatorial problem of searching among elements of the form 

f2le,Ie = ( U <5i) ^ ( U<5i), (10) 

and find the suitable index sets and Xq that minimize 

This simple idea is inspired by approximation theory of suitably expressing a func- 
tion as a linear combination of some given basis functions. Here instead of adding and 
subtracting the basis terms with suitable weights, suitable elements of 2) are combined 
through the union and relative complement to provide the shape approximation. 



We would like to highlight that more complex reference forms other than (10) may be 



considered. We however maintain simplicity by suggesting this form and assuming that 



D is rich enough that (10) still provides the desired flexibility in shape representation. 



For instance if the dictionary is poor, consisting of only two shapes Si and ^2, and the 



true shape is = 5i n 52, none of the possibilities in the form of (10) would be able to 
express (l. However if (l is already among the elements of the dictionary or we consider 
a richer dictionary consisting of Si, Ss = Si \ S2 and ^4 = ^2 \ 5i, we can find Ct by 



exploring different possibilities of the form (10) since - {Si u ^2) \ (53 u ^4). 



Selecting appropriate elements for 2) may be based on the level of prior information 
about the geometric features of fi. The shapes Si may be simple geometries that combine 
to form a more complex structure (e.g. see Fig[T]^a)); a large collection of shape possibil- 
ities among which the true shape needs to be determined (Fig[T]^b)); or a combination 
of both cases. 

To bring this idea into application we need to develop a computationally tractable 
algorithm that performs the search among possible shape combinations and minimizes 
the functional in ([T]) using a sufficiently sparse set of elements in the dictionary. More 
formally we reformulate ([T]) as 

{X@,Xq}= argmin "EiVtx^^Xe), (11) 

P®|+|Xe|<5 



where s represents the desired level of shape sparsity in the reconstructions. To pro- 
vide an approximate solution to the combinatorial minimization problem (11), in the 
sequel we present a relaxation strategy to maintain sufficient smoothness of the cost for 





(a) 



(b) 



Figure 1: (a) Reconstruction of an "Omega" -shaped region by applying basic set oper- 
ations on simpler geometries (two circles, a rectangle and a triangle): in this case the 
desired shape can be written as = {Si u ^2) \ (53 u ^4) (b) Choosing the right compo- 
nents of the desired shape among a dictionary of candidates: in this case extracting the 
word "SHAPE" from a collection of randomly placed characters 



gradient type optimization techniques to become applicable. We then provide a sparse 
reconstruction algorithm associated with nonlinear costs to employ a sparse number of 
dictionary elements in the final shape representation. 

2.1 Pseudo-Logical Shape Interaction 

In a recent work, Aghasi et al. provided a rather general view of parametric level set 
methods for inverse problems. As an example of a parametric form, they proposed a basis 
expansion of smooth compactly supported radial basis functions (named as bumps) with 
adaptive dilation and center points. They brought into attention a so called "pseudo- 
logical" property of this class of functions that we generalize its notion to broader appli- 
cability and something not specific to smooth radial bumps. 

Consider a given shape 5 c L) c R^. We use the terminology ^'knolF for a Lipschitz 
continuous function : ^ [0, 00) such that 

tps{x) > X G int(5) 
ips{x) = X e D \ S 

with int(5) denoting the interior of S. Intuitively, a knoll tps takes positive values inside 
S and vanishes outside S (and on the boundaries as a result of continuity). It can be 
easily inferred that for two knolls i/^Si and '0^2 corresponding to shapes Si and ^2, the 



(12) 



following facts are true: 

supp(^5i + '^^2 ) = <5i u 52 (13) 

and 

lim ^(supp^(^5i - c^^<S2),5i \ = 0. (14) 

Here supp"^(.) denotes the positive support, where the function takes values greater than 
zero and ^(., .) is a measure of dissimilarity between two shapes which basically vanishes 
for identical shapes (see for examples of this measure). The basic message is summing 
up two knolls would imply a union operation on their supports, and subtracting a knoll 
of large weight a » 1 from another knoll approximates applying relative complement on 
their positive supports. 

The idea may be easily incorporated with the notation of level sets by momentarily 
employing a c> level set instead of the zero level set. This unusual lifting is due to the 
compact support of the knolls, as using a zero level set causes ambiguity in identifying 
the underlying shape boundaries. In this context, for small values c > 0, the c-level set 
of i\)Sx +'0<S2 approximately represents S\ u52 and for large weights a » 1, the c-level set 
of '05^ - ail)s<2 approximately represents S\ \ ^2 (see Fig [2]). Reverting to (10), a more 
general corollary states that for a parametric shape defined as 




(a) (b) 

Figure 2: The pseudo-logical behavior of two knolls ^ind by considering a close 
to zero level set c (a) The c-level set of i\)s-y + '0^2 approximately represents S\ u ^2 (b) 
When ay>\ the c-level set of ^\ - a'02 approximately represents S\ \ ^2 



(15) 



making 

cjai 0"^ and ajjai +00, (for all i e e X©), (16) 

results 

^(f^?e,Xe,^^Xe,Xe)-0. (17) 

In other words, by considering the c-level set, the parametric function E^J^i is 
capable of producing shapes arbitrarily close to rix®,x©- To conform with prevalent level 
set methods we henceforth consider the zero level set by transferring the lifting parameter 
c into the level set formulation as 



(/){x,cx) = -c+Y^ai^/js,{x). (18) 
=1 



Of course combining the shapes by varying the coefficients a in (18) brings continuum 
into the problem. In other words, instead of solving a combinatorial problem and search- 
ing among a certain number of possibilities, our search is performed in a subspace spanned 
by '05- (x), which contains elements that can represent shapes arbitrarily close to rix®,x©- 
Searching in this subspace not only relaxes the combinatorial problem and makes it ca- 
pable of applying descent search methods, also provides the option of exploring shapes 
that are not among f^x®,Xe possibilities (Fig [3]). 

2.2 Generating Knolls from Given Shapes 

Earlier in Section 1.1, we pointed out the use of signed distance functions (SDF) in re- 
initialization of traditional level sets. In this context at every time step of the Hamilton- 
Jacobi equation the level set function is re-initialized as an SDF of the form 

^^^^^ [ -d(x,C) xeD^S ' 

where S is the underlying shape and d(x,C) is the distance between the point and the 
shape boundaries C. A similar idea may be used to uniquely generate knolls from a given 
shape. More specifically a knoll may be formed as 

, . [ d(x,C) xeS . . 

^^(^^ = [ xeD.S^ ^^^^ 

which basically implies '^^(x) = ps{x)^. Besides their widespread use in the level set 
community, our intention of using distance functions is to provide a rather basic ex- 
pression for the knolls and benefit the fact that shapes at the c-level sets of such knolls 
uniformly inherit many geometrical features from the original shape (see Fig [4]). More- 
over, since the the c-level set of a knoll il^s (or equivalently the zero level set of i/^s - c) 



is slightly smaller than 5, in a level set representation as (18), knolls may be generated 



with slightly larger supports, e.g., ips(x) = (^ps(x) + c) . By this choice, reconstructing 
exact elements of the dictionary does not require having very large weights associated 
with their knolls. 




Figure 3: Interaction of four identical 
circular knolls with different signs which 
gives rise to a triangular region with al- 
most straight sides. The resulting shape 
is certainly not in f^Xe,Xe possibilities 



Figure 4: A knoll corresponding to a star- 
shaped region. The shape resulted at 
some c-level set near zero inherits many 
geometrical features of the knoll's support 
such a the general structure and corners 



3 Sparse Reconstruction 



We so far discussed how considering a parametric level set as (18) relaxes (11) into a 
minimization of the form 

q: = argmin£^(a), (21) 

||a||o<s 

where S{cx) = ^{Q^f^^ f^^). We are certainly interested in sparse solutions of a to avoid 
shape redundancy and ill-conditioning associated with large dictionaries. In general £ 
is a nonlinear and non-convex function of a, we however assume it to be a sufficiently 
smooth function of ex. 

In the sequel we provide a general overview of sparse recovery techniques in find- 
ing solutions of underdetermined linear systems. We then use this notion to develop 
our sparsity promoting algorithm applicable to a large class of functionals appeared in 
imaging applications. 



3.1 Background on Sparse Recovery Techniques 



Finding sparse solutions of linear systems is a broad area of research in imaging science 
For an underdetermined linear system, Aot 



6, the main underlying problem 

(22) 



[28133143 

is 



mmimize 



s.t. 



Acx = b. 



The matrix A is m-hj-n where m « n and 6 is a vector of length m. Problem p2| ) is 
in general a hard combinatorial problem. However, it was brought into attention that 
relaxing (22) by replacing ||a||o with ||a||i, known as the basis pursuit (BP) problem 
can still result in sparse solutions. Moreover, under certain conditions on A, a 



BP solution perfectly coincides with the solution of (22) 12,13,22]. In case of noisy 



observations, the linear equality in (22) is replaced with a least squares inequality as 



mmimize 



ll<^lli 



s.t. 



\Acx - b\\2 < cr^ 



(23) 



where a is the noise level. This convex minimization is known as the basis pursuit 
denoising (BPDN) problem. We consider the BP problem as a specific case of BPDN 
when cr = 0. 

Another convex problem earlier used to induce sparsity on the solutions of a linear 
system is the least absolute shrinkage and selection operator (Lasso) ^45j, phrased as 



mmimize 



|Aa-6||2 s.t. 



< T. 



(24) 



For a certain value r = r^^ problem (24) becomes equivalent to (23), however, in general 



Tcr cannot be determined a priori. To promote sparsity, in most applications solving (23) 



is more desirable than solving (24). This is mainly due to the fact that a good estimate 



of the noise in the measurements is more likely to be available than prior knowledge of 
T, the ^i-magnitude of a sparse solution. 

A variety of solution strategies may be taken for both the BPDN and Lasso problems 
(e.g., see |11, 26, 34, 35, 45]). Relatively speaking, the smooth cost associated with the 



Lasso makes it an easier problem to approach compared to the BPDN which involves a 
non-smooth cost minimization. Minimization of smooth costs over convex sets may be 
efficiently handled using methods such as spectral projected gradient (SPG) |7 . 

We briefly overview the SPG-based technique proposed by van den Berg and Fried- 
lander in (47], who show that a numerically efficient way of solving the BPDN problem 
is by solving a series of Lasso problems. This technique plays a major role in developing 
the algorithm that we put forth later in this paper. Their method also allows use of A 
as a linear operator with no explicit matrix form which suits our general representation. 

Following [47], for a given r > 0, a single parameter function ^{r) is defined as 



(f{r) = \\b-A(Xrh, 



(25) 



where a^- is a solution of (24) for the given r. It should be noted that multiple solutions 
to the Lasso may exist, however, problem's convexity requires all of them to generate an 
equal error term (f^r). It is shown in |47| that is a nonincreasing convex function 
which is continuously differentiable and 

^'(r) = ^^ = -^^fc^^. (26) 
dr II 5 -Aa^ II 2 

The graph of ^(r) in terms of r provides an optimal trade-off between the residual error 
||6- Aar||2 and the ^i-norm of the solution, forming the Pareto curve (Fig [5]). The Lasso 
and BPDN problems are basically two different characterizations of this curve. In other 
words, for a given a a solution to 

^(r) = <7, (27) 
returns the value of Tcr that makes solutions of BPDN and Lasso problems coincide. As 



a matter of fact solving (27) does not require full access to the Pareto curve. The convex 
nature of the problem and a closed form expression for ^'{r) are sufficient to employ 
root finding techniques and acquire Tcr. In this context, a Newton iterative process 
^i+i = + At^5 where 

Ar, = ^^, £ = 0,1,.. (28) 

is capable of generating a sequence of Lasso parameters that superlinearly converge 
to Tcr. This result makes Lasso the central tool in solving the BPDN problem. More 
specifically, for a given a > 0, a solution to the BPDN problem is acquired by solving a 
series of Lasso problems parameterized by and ultimately arriving at a Lasso that is 
parameterized by Tcr, which basically solves the BPDN problem. To solve the underlying 
Lasso problems efficiently authors in |47] employ the SPG technique detailed in |7 . 

To address a broader class of problems, specifically in a variational framework, in 
our future formulations we look beyond matrix representatives of A. We consider a 
more general case that A : ^ S is a bounded linear operator and S represents a real 
Hilbert spac^ An example of S, other than R^, is the space of bounded continuous 
functions C(R). This generalization only requires slight modifications to the above 
overview of sparse recovery techniques, basically by replacing ||.||2 with ||.||§, the norm 
induced by the Hilbert space S, and using the adjoint operator in place of . The 
sparsity constraints applied to ol remain intact as the domain of A is still R^. For 
linear systems, the technique proposed in |47| makes such generalization possible and 
conveniently handles cases that A is not explicitly available as a matrix. 



^We therefore keep the notation general and do not use a bold syntax (representative of a matrix 
form) for the linear operator and elements of S. 



Figure 5: A typical Pareto curve corresponding to the Lasso and BPDN problems. For 
a given a, the corresponding value Tct makes Lasso and BPDN problems share solution. 
For values of r larger than tbp (the £i-norm of the BP solution), the residual \\Aa - b\\2 
is zero. 



3.2 A More General Sparsity Constrained Optimization 



The majority of sparsity-constrained optimization techniques are merely applicable to 
linear systems and corresponding quadratic costs. We are however interested in sparse 



stationary points of S{cx) in (21), which is not generally a quadratic function of a. 



This is a more complex problem and a new area of research which still requires further 



development. The current available techniques are mainly in the form of (21) for which 



we are required to know the degree of sparsity, 5, prior to the minimization [3|[4|| 

For the purpose of this paper, we consider a particular form of £{cx) applicable to 
a large class of imaging problems. We then provide a sparsity promoting minimization 
scheme that takes a rather short iterative process to converge and does not require prior 
assignment of s. Accordingly, for a given real Hilbert space S, we consider £{ot) to take 
a form (or a finite sum of forms) as 



S{a) = \\g{a)\\l = {g{a)Mc^)), 



(29) 



where ^ : ^ S is a Frechet differentiable map applied to a. The classic compressed 
sensing problem can be interpreted as a special case of (29) when Q{.) is an affine function 
of a, i.e., G{ct) = Aa. - b. In the following we present two classes of imaging applications 
that follow the form (29) and will be later considered in Section [4j 



3.2.1 Shape-Based Inverse Problems 



For this class of problems, in the basic case, the spatial parameter of interest u is modelled 

as 

u{x) = I ""^^ ' ^ ^ , (30) 

^ ^ [ Uex X e D \ U ^ ^ 

where Uin corresponds to the texture value of an inclusion (the shape) and represents 
the background texture. The texture values may be represented by scalars or low order 
parametric models which may/may not be known a priori 1 1 . Using a parametric function 
to characterize the shape by its zero level set we may rewrite the parameter of 
interest as 

u{x, a) = UinH[(j){x, ol)) + i/ex(l " H[(j){x, a))^, (31) 

where H{.) is the heaviside function, usually replaced with a smooth approximation 
Hrg{.) to maintain differentiability |1 . The inverse problem may now be formulated as 
minimizing the model-data mismatch to acquire the shape parameters a, i.e., 

mm\\v - M(u(x,ot))\\i . (32) 

Q, II \ / W^y 

For a smooth physical model A4 that relates the parameter of interest u to the observa- 



tions, the underlying inversion cost may be interpreted as a functional of the form (29) 
when G{ct) - v - A^(i/(x, a)). 

3.2.2 Image Segmentation 

A well known variation form to segment an image u{x) into two disjoint regions and 
D \ Q corresponds to a functional 

£(f^) = £,(J^)+ f r,n{u{x))dx+ f rex{u{x))dx, (33) 

where ^(fl) represents some regularity constraints on the segmented regions (e.g., 
smoothness and compactness |14, 16 ) and rin{.) > and rex{-) > are some inho- 
mogeneity measures of u in each region. A well known instance of such measures is 



observed in the work by Chan and Vese 16 , which suggests rin{u) - {u{x) - Uin)^ and 
T^ex{u) = {u{x) - Uex)^ ' Here Uin and 

Uqx are scalars representing the average pixel value 
within each region. Another example is using a maximum likelihood approach for the 
pixel intensities to arrive at 



{u) = -log(^p^{u{x))y X'in^ex, (34) 



where Pin{u) and Pexi"^) 9^re pixel intensity distributions [20,41 . The scalar parameters 



'Sm, Uex and the distributions Pin{-)i Pex{-) may be known a priori or may be estimated in 



the course of segmentation. Considering the image dependent terms in (33), a parametric 



level set function may be employed to form the segmentation functional 

£{cx) = rin{u{x))H{(l){x, a)) dx 

^ f^^ex{u{x))[l-H{^{x,cx))) dx. (35) 

We note that in using a parametric level set technique the regularity constraints may 
either be neglected thanks to the intrinsic smoothness of the parametric function [27] , 
or in our case taken into consideration via the knolls sparsity. The energy functional in 



( |35| ) may be written as sum of two functionals £in(ct) and Sexist) in form of (29) using 
Ginicx) = yJrin{u)H{(j){x,cx)) and Gexio^) = yjrex{u){l- H{(l){x,cx))). Employing a 
smooth function Hrg{.) : R ^ [0, 1] to approximate the Heaviside function can guarantee 
smoothness of Gin{-) and Qex{-) in oc. 

3.3 A Sparse Nonlinear Minimization Technique 

To minimize a sufficiently smooth cost function £{ol) : ^ R, a well known iterative 
scheme is Newton's method. Starting from an initial vector ag, Newton's method pro- 
ceeds by generating otk vectors {k - 1,2,-") that progressively reduce the cost to reach 
a minima. Having ock available, the successive vector is written as ock+i = cxk Sk (or 
a multiple of S^), where the step is determined by minimizing the second order Taylor 
approximation of the cost around a/^: 

Sk = argmin £{cxk) + Js{oLk)S + \6^1-Ls{oLk)d. (36) 

Here J's{oLk) denotes the Jacobian vector at OLk and Hs{(Xk) is the corresponding Hes- 
sian matrix. 



A positive definite Hessian guarantees existence of a minima for (36) which leads 
to the closed form expression Sk = {Jis{oLk))~^ J's{oLk) . Most Newton type techniques 
either impose such positivity constraint on the Hessian or approximate it with an at least 
positive semi-definite matrix to guarantee convexity of the underlying cost |21|. 

To consider applying a Newton type technique to the energy functional in (29), using 
calculus of operators [ll[29), we have 



Je{cx)5 = 2{g'{cx)5,g{oc))^, (37) 

and 

V^f(a)5 = (g'(a)5,a'(a)5)g + (g"(a)[^,^],^(«))s- (38) 

Here • R^ ^ S is a linear operator representing the first order Frechet derivative 

of ^ at a and Q''{ol)[.^ .] : R^ x R^ ^ § is a bilinear operator that represents the second 



order Frechet derivative of Q. When S is taken to be R^, and Q^^ are simply the 
Jacobian matrix and the Hessian tensor of G- For a positive semi-definite approximation 
to Hg^ in Gauss-Newton methods |21| the inner product term containing Q^^ is neglected 



in (38). Therefore the underlying quadratic cost corresponding to (|36[) becomes 



\\g{cxk)\\l + 2 {g'{ak)s, QM)^ + \\Q'Ms\\l 

= \\g'{ak)S + g{ak)\\l 

forming a convex function of S. Although the underlying quadratic cost is convex, the 
linear operator may be ill conditioned and some type of regularization may be required 
to obtain S^- For instance Levenberg-Marquardt algorithm is a variant of Gauss-Newton 
technique that uses a Tikhonov type regularizer to minimize the underlying quadratic 
cost |21|. 

We however consider a different regularizer which is capable of promoting sparsity 
on consecutive estimates of ock- More specifically having otk available we regularize the 
subproblem by promoting sparsity on the next potential estimate cxk+i through a BPDN 
problem, i.e., 

Sk = SiTgmm \\cxk-^S\\i s.t. \\G\otk)S -\- Q{cxk)h < (39) 
d 

This is certainly a translated version of the BPDN problem and a simple change of 



variable as = a/^ + 5 converts it to the standard form ( |23[ ). As S is expected to be very 
small about an accumulation point, a reasonable choice of a is an estimate of ||^(oi)||g 
near a minima. In Appendix A we show that when a = min^, ||^(a)||g, or even more 
locally when a < \\Q{cXk)\\'^, the direction acquired in (39) is a descent direction. 



In the case of inverse problems, a may be an estimate of noise power in the measure- 



ments. Nevertheless, if such estimate is not available, a substitute for (39) would be a 
parameter free BP problem 

dk = argmin^ ||ajt + 

s.t. g''{cxk)G\cxk)s + g'^MOM = o, (40) 

corresponding to the classic Gauss-Newton equation with an additional sparsity con- 
straint. If Q\cXk)S + G{otk) = has a class of solutions, which is very likely for an 
ill-conditioned linear system, this equality constraint may be replaced with the one in 



(40) and make it equivalent to the case a = in (39). 

As stated in Section 3.1, at every iteration a minimizer to (39) may be determined 
via a sequence of Lasso solves 

Sk/ = argmin^ \\g\(Xk)S ^ G{ock)\\s 

s.t. \\cxk + S\\i<Tk/, (41) 



where for i = 0, 1,2,--- the ii radius r^/ is updated according to (28) and sequentiahy 
hm^^oo = Sk- Once the fc-th iteration is complete, a similar Lasso search is performed 
along the Pareto curve corresponding to \\Q\cXk+i)S + G{cxk+i)\\s to find S^+i- 

In determining the consecutive directions dk, it is computationally desirable to per- 
form incomplete iterations on £ and make a more accurate search for when a.j^ ap- 
proaches a sparse solution. For this purpose, it is also shown in Appendix A that if for 
some ^, Tk/ > \\a.k\\i, the corresponding step, Sk/^ is a descent direction and can be used 
as a step to reduce the nonlinear cost. 

Also, for small values of k that a.k may be far from a sparse solution, Tk/ may tend 
to very large values because of the high cardinality of otk- To prevent this, the values Tk/ 
may be controlled by a loose threshold value Tmx and the iterations on £ may be broken 
if this threshold is reached. The descent property of the incomplete step for this case 
is also inferred by a similar argument as the one just mentioned. These two strategies 
provide us with an indication of when we are eligible to break the ^-iterations and still 
be certain that the descent property of resulting direction is maintained. 

Algorithm 1 provides a detailed picture of the proposed minimization scheme where 
a line search technique is employed to warrant a proper convergence behavior. We 
specifically use Armijo rule to determine the step size, which requires the reduction in 
the cost to be sufficiently large |6 . Parameters 7 and /3 correspond to this reduction 
and usually take values as 7 g [10"^, 0.1] and /3 g [0.1,0.5], detailed in [g]. The values 
£1 and £2 control the size of updates in each block and are very small numbers in the 
order of machine precision. The value of Tmx be taken as a loose overestimate of 
the ^i-norm of the solution. Finally cxq is the initialization of unknown parameters and 
To is an initial ^i-ball radius which may arbitrarily taken to be zero or very small. 

3.4 Employing an Asymmetric Norm 

In the sparsity promoting algorithm proposed we used the £1 norm as a convex relaxation 
to the £q norm. This is a common technique in compressed sensing promoting sparsity 
on the reconstructions thanks to the sharp vertices of the £1 ball. However, in general 
an £q minimization only targets the cardinality of a vector while an £1 minimization also 
takes into account the vector component values. 



For a parametric representation as (18), Fig[6[a) shows reconstruction of an L-shaped 
region where a basic £1 minimization does not necessarily provide the best sparse solution. 
Considering a dictionary that contains all shapes Si,"',Se shown in the figure, in one 
case the L-shape is reconstructed via the union of four rectangular knolls while a similar 
reconstruction is possible by applying the relative complement among only two knolls. 
When ||a||i, the £1 norm of knoll coefficients, is used to indicate the level of sparsity, 
the latter case may have a lager £1 norm although only two shapes are exploited. This 
is because of the relatively large negative coefficients that are required to simulate a 



Algorithm 1 A sparsity promoting Gauss-Newton algorithm 
input cr; 

set 7, /3, 51, 82 and Tmx] 
OL := ao; 

T := To; 

found := false; 
while ~ found do 

step found := false; 

while ~ step found do 

d:=argmin^ ||^^(a)r7 + ^(a)||§ s.t. ||a + ?7||i<T; 
<^;:=-||a'*(a)[e'(a)5 + G{a)] \\j\\G'{c^)S + 
At := {a-\\g'ia)5 + g{a)\\s)/<t>U 
r := r -\- At; (alternatively, r := min(T + At, Tmx)]) 
step found := (t > ||a||i) v (At < £i)] 
end while 

while £(a.) - £{a, -\- S) < -^Js{ol)5 do 

5:=/3>^5; 
end while 
a := a + (5; 
found - (||(5||2 <£2); 
end while 




(a) (b) 

Figure 6: Employing an asymmetric £i norm: (a) Reconstructing an L-shaped region by 
union of four rectangles (left) and applying the relative complement among two rectangles 
(right); (b) The £i ball and the asymmetric convex ball corresponding to ||oi|i,i(; 

relative complement operation (revisit Fig|2]). 

One way to remedy such phenomenon is to consider an asymmetric £i norm defined 

as 

ll<^li,^= X! Z! M^^jl (42) 

{i:ai>0} {j-Oij<0} 

where it; is a constant scalar. When it; g (0, 1), this representation allows the appearance 
of larger negative coefficients by weighting them less in the absolute sum. Addition- 
ally, setting w » 1 promotes positivity on the coefficients which is desirable in some 
applications that we will consider in next section. 

Clearly \\.\i^w is not formally a norm as it violates the positive homogeneity of a norm, 
however, the corresponding ball is still convex (see Fig[6]^b)). A generalization of |47j 
reported in |48 allows replacing the £i constraint with any convex constraints. As stated 



in [48] for this generalization the infinity norm in (26) needs to be replaced with its polar 
(reducing to the dual norm when the constraint is a norm), which in this case is the 
asymmetric infinity norm 

||a|oo^-i= mdix {\ai\,w~^\aj\}. (43) 

' i: ai>0 



To more conveniently facilitate our method with this feature, we would note that in a 



lower implementation level, a main component of the SPG method in ^47j is performing 
iterative projections of the form 

a"^ = argmin ||a-a||§ s.t. \\ck.\i^^<r. (44) 

The asymmetric norm maybe written as ||o:|i^^(; = ||Z)(a)a||i where D{cx) is a diagonal 
matrix with entries 

^M,. = |' ""'^n- (45) 



As we have shown in Appendix B, solution of (44) coincides with the solution of the 
weighted £i minimization problem 

min ||a-a||§ s.t. \\D{a)cx\\i < r. (46) 

This fact is well demonstrated in Fig [7| Based on this argument we may replace the 
projection onto an asymmetric ball problem with a weighted £i minimization, which is 
straightforward. 




Figure 7: Coincidence of the projections onto the asymmetric balls ||oi|i^-i^ < r and 
||£)(a)a||i < T. The two balls share side in the second quadrant where a is located. 



4 Examples 



In this section we examine the proposed technique in various imaging applications. For 
the examples presented, we use the approximate Heaviside function suggested in |1| as 



H,{x) 



1 




X > e 
X < e 
\x\ < e. 



(47) 



which gives rise to a compactly supported approximation of Se{x) = H'^{x). As discussed 
in [l] in the context of inverse problems, this choice results in a form of narrow-banding 
in the evolution of the parameters, that will also be discussed briefly in segmentation 
examples. The parameter e determines the width of the narrow-band region, whereas 
smaller values of e brings sharper transitions into the reconstructions at the expense 
of slower evolution |1|. The choice of the lifting parameter c is quite arbitrary as the 
proposed shape-based problem is relative, i.e., simultaneously scaling c and ai coefficients 



in (18) does not change the zero level set shape. In other words, taking larger values of 
c would tend to larger values of ai in the reconstructions. Since in general the shapes 
in the dictionary may significantly vary in size, to have knoll basis terms of comparable 
magnitudes, we normalize each knoll to its maximum value. For the majority of examples 
we simply take c = 0.1 and e = 0.05 and initialize the algorithm from a rather random 
state. Based on the scaling property between c and final ai values, r^x i^iay be taken as 
a large multiple of c, e.g., Tjjix ~ 

50c. 

4.1 Image Segmentation 

As discussed in Section [3.2.2| for a binary Chan-Vese segmentation the energy functional 



finds the form of (35) which maybe written as the sum of internal and external energies: 

S{oc) = \\Qin{oc)\l+\\Qe:c{cc)\l (48) 

Accordingly a reasonable representation for Q{ol) would be 

G{0c)=G^n{0i)+xGex{0c), (49) 



where i = The underlying Hilbert space would be § = L?{D) with the corresponding 

inner product 

52(x))^2(2^) = 5l(x)52(x)dx, (50) 



where 52 (x) represents the complex conjugate of 52 (x). Clearly, applying our proposed 
algorithm to the segmentation problem requires having ^^(o:)[.] and its adjoint form. 



For a given vector 77 g M^^, a formal derivation of (49) with respect to ol and applying 
the resulting Jacobian operator to ry yields 



rid 



g'{oL)r] = /(x,a) ^77i^5,(^), (51) 

i 

where 

/(.,c) = | ^(V^-'VSij)*'**' 1*1^^ (52) 

else. 



From a numerical perspective, the linear operator G\cx) is a matrix with columns 
f{x,cx)t(jSi{x). Based on this argument, the adjoint operator may be easily deduced. 
For a function e L^(D) the linear operation Q'''{cx)i/j produces an array of size 
where 

[g'*{oc)i;]i = (J{^i;s.{x),i;{x))^^^^^. (53) 



To maintain generality, we preferred providing the implicit forms (51) and (53) for G\cx) 
and ^^"^(a), however, for pixelated images of moderate size a matrix may be assigned to 
each linear operator. 

Equation (52) shows that /(x,a) is a compactly supported function. Similarly, 
'05. (x) is compactly supported, and therefore when supp(5e((/))) n supp('^/^5. ) = 0, the 
corresponding columns of G'{cx) would become zero. In determining a descent direction, 
such columns may be neglected as they would not affect the result. In other words, in it- 
eratively updating the knoll coefficients a^, the only knolls that are updated are the ones 
that intersect the narrow-band 5(0) at that iteration. This numerical advantage is basi- 
cally the narrow-banding feature associated with the proposed parametric representation 
and extensively described in |1 in the context of inverse problems. 

The technical details presented may be employed to address variants of Chan-Vese 
segmentation in different applications. Here we consider two applications that could be 
addressed through a sparse shape recovery. 



4.1.1 Image Segmentation with Missing Pixels 

A challenging problem associated with image segmentation arises when some of the pixels 
are missing or part of the image is occluded. This limits the analysis to only a subset 
of image pixels, while the segmentation needs to be performed globally. In this example 
we show that prior information about the geometry of objects in the image can lead to 
perform a completion to missing pixels when a sparse shape recovery is considered. 

Fig |8] shows a noisy reference image that is occluded with rectangular patches. Ad- 
ditionally, random pixels of the image are discarded resulting in overall pixel loss of 50% 
and 80% shown in Figures [sj^b) andjsj^c). These two figures will be the subject of our 
segmentation. 

Our prior information about the geometry of objects in the image is reflected in the 
choice of dictionary elements. To build up the dictionary we make use of four basic 
shapes: a circle, square, triangle and an ellipse, shown in Figjsj^d). We would note that 
the triangle in the reference image is the upper diagonal portion of a square while the 
one used in the dictionary is a lower portion. 

A slight modification of the formulation presented above needs to be considered for 
this limited observation problem. More specifically, considering c D to be the available 
portion of the imaging domain the Hilbert space considered would be L^{D') where 





(d) (e) (f) 

Figure 8: (a) Reference image (b) A test image with more than 50% of the pixels missing 
(blue regions show the missing regions) (c) A test image with more than 80% of the pixels 
missing (d) Shapes used to build up the dictionary: 1000 instances of these four shapes 
with different sizes are placed throughout the imaging domain (e) Initial values of ai for 
i = 1,-"1000 (f) The resulting initial contour for the given initial ex. The segmentation 
results are shown in Fig [9j 




Figure 9: First row shows the results of segmentation on 50% missing-pixel image, from 
left to right for it; = 0.1, it; = 1 and w = 10. Below every result the final sparse vector a is 
shown. Second row are the level set functions corresponding to each segmentation above 
it. Third row shows an identical scenario as the first row applied to the image with 80% 
missing-pixel 



the observation takes place. The knoh functions, however, are constructed and defined 
globally in D. 

By considering different size and placements of the four basic shapes, a total of 
rid - 1000 shapes are assigned to the dictionary. To initialize the algorithm we randomly 
apply weights of +1 and -1 to 100 knolls as shown in Figjsj^e). This assignment results 
the initial shape shown in Figjsj^f). A basic Chan-Vese segmentation cost is considered 
for which the values uin and Uext are occasionally updated following the formulation 
in ^6|. 

The first column of Fig |9] shows the successful segmentation results on the two test 
images that are obtained for it; = 0.1. Although significant parts of the images are 
missing, in both cases the proposed algorithm has made a reasonable segmentation job. 
To highlight the advantage of considering an asymmetric norm, results are also shown 
for the basic case of it; = 1 and an inappropriate choice of it; = 10 which basically 
abandons set minus operation. It can be seen that although considering smaller values 
for w increases the size of feasible region and requires solving the problem in a larger 
domain, it pays off by promoting the set minus operation and using a sparser set of 
shapes in a more efficient manner. This contrast is clearly observable by comparing 
the results of first and third columns. Basically, in the third column knolls are pushed 
to take positive weights which degrades the reconstruction by exclusively promoting the 
union operation. Thanks to the appearance of the set minus operation, the segmentation 
results in the first column have sharper corners and smoother sides closer to the truth 
while using relatively less number of shapes in the dictionary. It is worth noting that all 
segmentations converged the steady state in less than 15 iterations. 

4.1.2 Text Recognition: Breaking a Basic Captcha 

Captcha (Completely Automated Public Turing test to tell Computers and Humans 
Apart) is a well-known test in the world of computers that ensures that response to 
a query is generated by a human |50 . For this test, the client is asked to read and 
type the word in an image, where the characters are placed in an unusual manner not 
recognizable for the machine. Since the underlying components of the image are still 
characters, this maybe considered as a challenging sparse shape reconstruction problem. 
The underlying shapes in the dictionary may constitute a dense set of possibilities for 
the character shapes that may appear in a Captcha image. Of course a denser dictionary 
increases the chances of identifying the characters correctly. 

As a proof of concept, in this example we consider an attack on the rather basic 
Captcha image of size 63 x 160 pixels shown in Fig [To|(a) which represents the word 
"S'/iaPE"'. The letters in the image may take different case, size, rotations and overlap as 
the case for our Captcha image. To build up the character dictionary we assume knowing 
the font type a priori and consider all 52 uppercase and lowercase letters of the English 



alphabet. The displacement possibilities considered for each character are 64 points on a 
regular grid in the image. At every point, five different rotations of the character and two 
different font sizes are considered. This setup leads to = 52x64x5x2 = 33280 shapes in 
the dictionary. Of course for a more realistic scenario, font possibilities, deformations and 
more number of size, displacement and rotational variations may be considered for every 
letter. We however keep the problem small to be tractable with a desktop computer. 
Since the main purpose of this problem is identifying shapes in the dictionary that appear 
somewhere in the image, w may be chosen to be large (10 for our simulations) to push 
the algorithm on only considering the union operation. 

As the first experiment we consider the case that exact character shapes in the 
Captcha are present in the dictionary. As before, the algorithm is initialized with pos- 
itive and negative weights for a random subset of ai coefficients which corresponds to 
the initial shape contour in Fig 10 ^b). The classic Chan-Vese cost functional is again 
considered with the texture parameters updated occasionally. The segmentation result 



after 12 iterations is shown in Fig[TO[c) and in Fig 10 ^d) we have shown the coefficients 
values. A simple index map relating the weight indices to the alphabetical representation 
of the shape indicates that the reconstructed image is composed of five main characters, 
precisely matching with the word inside the image. 

For the second example we consider the Captcha in Fig [TT](a) where the letters in 
the image do not precisely match the dictionary elements. The closest elements of the 
dictionary are depicted with colored contours on the same image. This problem may be 
still solved using a multi-stage refinement strategy. More specifically, a first segmentation 
attempt shown in Fig [TT](b) with the underlying weights shown in Fig pT|c) does not 
provide a sparse solution with prominent weights as before. Instead, in this case there 
are several candidates, the top 20 of which are shown in Fig[TT](c). Based on the first 
stage results, a new dictionary maybe built up, with much a much smaller number of 
elements, only containing denser size, rotation, displacement variants of the shapes listed 
in the top 20 list. Performing the segmentation process again, this time provides us with 
the weights shown in Fig pTj^e) where the striking weights only correspond to the letters 
in the Captcha. The "ec/io" effect corresponding to appearance of multiple prominent 
weights of the same letter is due to the close distance between the shapes in the new 
dictionary. Of course when the indexing mechanism takes into account the character 
positions, all such echoes correspond to close positions and maybe unified as a single 
final conclusion about the letter in that region of the image. 



4.2 Medical Imaging: X-ray Computed Tomography (CT) 

As a well known linear inverse problem, in this section we examine the method in recon- 
structing CT images. For a mono-energetic CT, X-ray photons are transmitted through 
the test medium and measured at the opposite side. If the medium has an attenuation 
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Figure 10: (a) A Captcha image (b) The initial contour corresponding to the random 
initiahzation of the ai coefficients (c) The segmentation result shown with colored con- 
tour (d) The values of reconstructed weights and an indication of the letter each index 
corresponds to 
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Figure 11: (a) A Captcha image where the letters in the image do not match with the 
elements of the dictionary; the closest elements of the dictionary are shown with colored 
contours (b) A first stage reconstruction results (c) Reconstructed weights corresponding 
to the first stage segmentation; the top 20 candidates are indicated with the correspond- 
ing letter (d) Segmentation results after refining the dictionary according to the first stage 
results (e) The weights corresponding to the second stage reconstruction indicating the 
letter each index corresponds to 



profile the number of photons measured would ideally be 



A. 



ATexp(- / ^(x)dx), m = l,---,M, (54) 



where is the photon count measured at m-th receiver, A^^ is the blank scan photon 
count and Cm is the line through which the ray travels. An easier way of interpreting 
the measurements is reading the values 

i;^ = -log^= [ ii{x)dx, (55) 

at each sensor which basically relates the attenuation map to the measurements via a 
Radon transform. In practice however, the measurements are corrupted with Poisson 
noise, i.e., 

A^ = Pois(A^) (56) 

and the true measurements are = -log(A^/AT). A second order approximation to 
the log likelihood function indicates that |44) 



\ogp{v\lJi) ^--{v- niifD{v - TZfi) + h{v) (57) 
2 

where D = diag(Ai, Am), is the linear Radon type transform that maps fi to the ideal 
measurements and h{.) is a function dependent only on the data vector v. To apply the 
proposed algorithm to this modality the attenuation profile is parameterized as /x(x,a), 
and the residual operator is written as G{ct) = 7^/x(x, a) - v. For a maximum likelihood 
estimate of the parameters, the underlying cost takes the form ||5(o^)|||) where the inner 
product in the measurement space is defined as (si, 82)0 = sJDs2 for si, S2 ^ M^. Based 
on the linearity of the CT model the Jacobian operator applied to a vector 77 may be 
written as 

g (a)r, -T.V^^- T^V^n— = (58) 

and the adjoint operator applied to a vector -0 g takes the form of 

g^^(cx)^|, = {n^,^l^)n. (59) 
dai 



Equations (58) and (59) are in general the key components of applying Algorithm 1 to 
this problem. 

Fig [12]^ a) shows the true attenuation profile of a 200 x 200 pixel chest test image. To 
generate the X-ray data the rays are emitted at 60 equispaced angles between and 180 




Figure 12: (a) The true attenuation profile of a chest image (b) Inversion using FBP 
(c) TV inversion (d) Sparse shape inversion results (e) The basic shapes used in the 
dictionaries 



degrees. At each angle an average of 38 photon rays are transmitted through the image. 
A total of 2280 measurements are obtained among which 1499 are the ones crossing the 
chest and containing useful information. The blank scan photon count is = 4 x 10^ 
and the average photon energies are 50 KeV. Beside considering the Poisson noise, 1% 
zero mean Gaussian noise is also added to the data to model the inaccuracy in the 
measurement readings. This problem setup poses a challenging ill-posed problem which 
is rather hard to approach. In Fig [T2[b) we have shown the filtered back projection 
(FBP) results for this data set. Fig |12[c) also shows the reconstruction results using a 
total variation (TV) approach employing the ^i-magic package [TT] . 

For the purpose of shape representation we consider a multi-phase level set approach 
[49] . More specifically to invert for both the soft tissue and bone geometries along with 
the blank space we consider using two level sets (/)i(x, a^-*^^) and (/)2(x,a(^^) as 

fi{x, a(i), a(2)) = i_i, + - ;x„)i?,(</.i(aW)) 

+ (/ife - Ai.)i^.(</>i(a«))if,(,/.2(a(2))) 

where /i^, /^s and /i^ are respectively the average attenuation values for air, soft tissue and 
bone. Considering the average density of each material, at 50 KeV photon energy, the 
average attenuation values for the air, soft tissue and bone are approximately 2.7 x 10~^, 
0.2 and 0.7 cm~^ |30j. We use two shape dictionaries 2)^^^ and The former contains 

n^J'^ = 5346 shapes and the latter is in hold of rSP = 6156 shapes. The basic elements 

Elements of D^^"^ are chosen to 



12 



used in the dictionaries are those shown in Fig 
have larger sized shapes as 0i is mainly in charge of representing the soft tissue. It 
is however worth noting that both dictionaries share identical shapes. The shapes are 
placed all over the imaging domain specially in places that there are chances of objects 
being present. Plugging the sensitivity relations 



and 



da) 



(2) 



into (58) and (59) provides the necessary components in using Algorithm 1 to invert 
for the level set coefficients. The process of determining a descent direction may be 
performed in a coordinate descent fashion by alternatively updating the weights for the 
first and second level set at each iteration. Fig 12 ^d) shows the result of our inversion. 
Through a sparse composition of shapes we have been able to reconstruct a reasonable 
estimate of true image. Considering the image to be a composition of shapes provides 
a better pose to the problem compared to the case of considering it as an image with 
piecewise constant regions suitable for a TV inversion. 



4.3 Resistance Tomography: An Archaeological Problem 



As the last example we consider an archaeological application of imaging subsurface 
tombs using electrical resistance tomography (ERT). The severely ill-posed nature of the 
problem as well as the limitations in placing the measuring sensors make this problem 
very challenging specially when the structures are closely spaced [25]. 

In this technique electric current is injected into the ground and some sensors measure 
the resulting potential on different regions of the imaging domain. Based on these po- 
tential measurements an inverse problem is solved to reconstruct the profile of electrical 
conductivity. The governing physics may be described as 

V • (crVp) = j in D, 

Ci^^+6p = ondD, (60) 

where_/ denotes the pattern of injected current, a is the conductivity and p is the resulting 
potential. Functions and ^2 are functions defined on the boundary of the imaging 
domain and are in charge of imposing appropriate boundary conditions [2]. 

To maintain simplicity, consider ps to be the potential resulted from a point source 
current j{x) - 5{x - Xg) and pg^m representing the measured voltage at points Xm in 
the domain for m = 1,---M. The forward model is a nonlinear operator that maps 
(j{x) to the potential measurements. To apply the proposed technique we need to know 
about the model Jacobian operator. Using the adjoint field technique, it can be shown 
that perturbations in the measurements are related to the conductivity perturbations 
through [2j 

5ps,m = SaVps ' ypmdx. (61) 

Here pm is the potential resulted from placing the point source current at x^, known as 
the adjoint field. When a parametric form a = a) is considered for the conductivity, 
we have Sa = X!i Saida/dai, where Sai is the perturbation of the z-th element of a. Based 
on this argument 

r 0(7 

Q'mi0i)V= j^^Ps-^Pm{Z^i-Q^)dx. (62) 

For a G M^^, the model Jacobian operator G\oc) is a matrix of size Mxri^ the (m, i) entry 
of which is jj^da jdaiV ps ' S/Pmdx. Accordingly, the adjoint operator Q^'^{cx) is simply 
the transposed matrix in this case. This notion can be extended to more complex cases 
of running different experiments with different sources, as well as more sophisticated 
source settings such as electric dipoles. 



Fig (13) shows the setup for a field experiment in a region where two closely placed 
tombs exist. A total of 50 sensors are placed on the ground and 20 experiments are 
carried out at each experiment using a pair of sensors as the electric dipole (shown with 



Figure 13: The ERT setting and sensor configuration 



identical numbers) and the remaining sensors as the measurement sensors. The data is 
polluted with 1% white noise. The true cavities are modeled as cubic structures shown 
in the figure. The structure on the left is slightly tilted in both azimuth and elevation. 



Following (31) the conductivity distribution is modeled as 

cr(x, a) = + (cTa - as)He{(l){x, a)), (63) 

where a a and ag are the average known values for the conductivity of air and soil. A 
total of rid = 14157 shape are used in the dictionary. The dictionary elements are cubic 
structures distributed all over noting that the true cavities are not among the elements 



of the dictionary. Fig 14 ^b) shows the initialization of the algorithm and Fig[T4[c) is the 
imaging results after 11 iterations. 

Reconstruction results after 43 iterations using conventional level set technique are 
also shown in Fig[l4|(e). We would like to note that the conventional level set inversion 
is initialized with a very good initialization as shown in Fig p^e). As observable, the 
proposed technique is able to provide a more accurate profile of the subsurface conduc- 
tivity and successfully separate the two structures. The counterpart, however, is unable 
to provide such level of detail and fails to make the separation in deeper regions of the 
ground that the sensitivity values are lower. 



5 Conclusion 



The idea presented in this paper may be considered as a new technique in approaching 
variety of shape-based imaging problems. The main message of this work is changing the 




Figure 14: (a) The true tombs (b) Initialization of the algorithm (c) Reconstruction re- 
sults (d) Initialization of traditional level sets (e) Reconstruction results using traditional 
level set technique 



geometric problem of shape composition into a variational problem that may be more 
conveniently analyzed. This conversion is performed rather simply through the notion 
of pseudo-logical property. In fact following this idea, the shape composition problem 
becomes similar to the classic problem of representing a function with a weighted sum 
of functional elements from a dictionary. The notion of sparsity comes next as a means 
of choosing proper elements, however, the shape-based nature of the problem requires 
struggling with a nonlinear problem. By considering several examples from different 
applications, we showed that, although a nonlinear problem, the proposed sparsity pro- 
moting technique can successfully handle struggling problems with large dictionary ele- 
ments. However, more general techniques applicable to larger class of functionals without 
restrictive assumptions such as knowing the sparsity degree a priori are desirable and 
welcomed. Of course the non-convex nature of the shape-based problem does not allow 
talking about uniqueness of reconstructions, however, our group is still working on ap- 
propriately restricting the problem for which such analysis is possible. This is still an 
open arena of research that requires exploring new techniques in nonlinear problems with 
sparsity constraints. 



A Descent Property of Proposed Steps 



We first show that when a is sufficiently small a direction acquired from (39) is descent. 
First a quadratic expansion yields 

\\G'{ock)5 + Qicxk) Hi = WQicxk) Hi + Je{a)S + \\g\cxk)S\\l 

As 6k needs to meet the constraint in ([39]) we must have \\Q\a.k)Sk-^G((y-k)\\E> ^ o-, which 
requires 

JeMSk <a- \\G{cxk)\l - \\g'{cxk)6k\l 

^^-mccu)\l (64) 



The most right side of (64) is certainly negative for sufficiently small values of a, when 
for instance a - miuQ, ||5(Q^)||g or more locally a < \\Q{a.k)\\'^. Such choices make Sk a 
descent direction for S Sit otk- 



We next show that when rk£ > \\cxk\\i, direction Sk£ acquired from (41) is descent. 



Since Ha^^^^ + 0||i < Tk/ and Sk/ is a minima for (41) we need to have 

\\G\ock)Sk/^GM\\l<\\G\ock)o^g{cxk)\\l 



II 2 
IIS' 



which using (64) simphfies to 

Js{ock)6k/ < -\\g'{ak)Sk4l < 0. 



(65) 
(66) 



B The Equivalent Problem to the Projection onto an Asym- 
metric £i-Ball 



Proposition: For a g given, consider aj; to be the projection onto an asymmetric 
£i-ball as 



a]^=argmin ||a-a||§ s.t. \\c^\i^w^^i 
and to be the solution to the weighted ii projection problem 
a2=argmin ||a-a||§ s.t. ||Z)(a)a||i < r, 



(67) 



where the diagonal matrix i?(a) is defined as (45). Then 0.1=0.2- 



Proof: We first note that both problems ( |67[ ) and ( |68| ) are framed as a projection 
onto a convex set, and in both cases the minima is unique (see Prop B.ll in |6j). Clearly 
elements of Oi and o must have identical signs over the support of a|; otherwise by 
changing the signs of incompatible elements a feasible point is obtained with a lower 
cost, which contradicts o^ being the global minima. A similar argument holds for 02 
and a, and so 

D{ai)oi = D{a)o^, (69) 



and 



D{a^)o^ = D{a)o^ 



(70) 



Suppose that 0.1^0.2. Since O2 is a solution to (68), it needs to satisfy ||Z)(a2)a2 ||i < r. 
This result along with (70) and the fact that ||q^|i,w = ||i?(a)a||i, reveal that O2 is also a 



feasible point for problem (67). Therefore, based on the strict convexity of the cost 



(71) 



In a similar fashion equation (69) reveals that a| is a feasible point for problem (68) and 
therefore ||a2 - 6:||§ < ||aj; - a||§ which contradicts (71). 
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